	PRO ICREC_GEM,Lx,Lz,lambda,apert,kx1,kz1
; Initial conditions for magnetic reconnection 2.5 D
; GEM Challenge. Double precision numbers
; Input: N = # points, lambda = width of the current sheet
; apert = amplitute of perturbation with kz1 and ky1 

N = 512
x = (Lx/511)*dindgen(N) - Lx/2d  ; periodic boundary conditions
z = (Lz/511)*dindgen(N) - Lz/2d  ; periodic boundary conditions
Nx = n_elements(x) & Nz = n_elements(z) & B0 = 1d
loadct,3 ; paleta de colores

bxdf = dblarr(Nx) & bx = dblarr(Nx,Nz) & a = dblarr(Nx,Nz);

for j = 0,n_elements(z)-1 do begin
for i = 0,n_elements(x)-1 do begin
	pert = apert*cos(kz1*z(j))*cos(kx1*x(i))
	bx(i,j) = B0*tanh(z(j)/lambda)
	a(i,j) = B0*lambda*alog(cosh(z(j)/lambda)/cosh(Lx/(2d*lambda))) + pert
endfor
endfor
print,bx[0,0],bx[0,n_elements(z)-1]

bxdf(1:N-2)=(a(0,2:Nz-1)-a(0,0:Nz-3))/(z(2:Nz-1)-z(0:Nz-3)) ; bx from fin. diff.
bxdf(0)=(a(0,1)-a(0,Nz-1))/(z(2)-z(0)) ; first one 
bxdf(N-1)=(a(0,0)-a(0,Nz-2))/(z(2)-z(0)); last one 

window,2,xsize=300,ysize=300
plot,z,bx[0,*],xr=[-Lx/2.,Lx/2.],yr=[-2.,2.],linestyle=1,xsty=1,ysty=1,xtit='y',ytit='Bx(x,y)',tit='Initital Conditions of B(x,y) using finite differences',col=black,back=1
oplot,z,bxdf

window,1,xsize=300,ysize=300
plot,z,a[0,*],xr=[-Lx/2.,Lx/2.],yr=[-20.,20.],xsty=1,ysty=1,xtit='y',ytit='a(x,z)',tit='Initital Conditions of a(x,z)',col=black,back=1

Nlev=10
alevp= 0 + (1+findgen(Nlev))*(max(a) - 0)/(1 + Nlev)
alevn=min(a)+(1+findgen(Nlev))*(0 - min(a))/(1 + Nlev)
window,3,xsize=Nx,ysize=Nz,title='a(x,z)'
contour,a,x,z,lev=alevp,xr=[-Lx/2d,Lx/2d],yr=[-Lz/2d,Lz/2d],xsty=1,ysty=1,xtit='x',ytit='z',tit='Initial Countour Labels of a(x,z)',col=black,back=1
contour,a,x,z,lev=alevn,c_linest=1,/overplot
contour,a,x,z,lev=0,th=2,/overplot

; Calculo y grafico: J(x,y)
jj=dblarr(Nx,Nz)
k2=dist(Nx)^2 
ak=fft(a,-1)
jj=float(fft(k2*ak,1))
window,4,xsize=N,ysize=N,tit='j(x,y)'
tvscl,jj

meanb = mean(bx)
meana = mean(a)
print,'Valor Medio a(x,y),B(x,y): ',meana,meanb
print,'a(x=0,y=0),a(x=0,y=2*pi) ',a(0,0),a(0,N-1)
print,'b(x=0,y=0),b(x=0,y=2*pi) ',bx(0,0),bx(0,N-1)
print,'bxdf(y=0),bxdf(y=2*pi) ',bxdf(0),bxdf(N-1)
end
